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The  flow  of  air  to  feed  oxygen  to  the  cathode  of  each  plate  in  a  proton  exchange  membrane  fuel  cell 
(PEMFC)  is  studied  for  a  300  W  stack  in  a  realistic  3D  configuration.  Two  configurations  for  gas  income 
are  solved,  a  “U”  shape,  where  both  the  inlet  and  outlet  of  the  air  collectors  are  at  the  same  end  plate,  and 
a  “Z”  shape,  where  inlet  and  outlet  are  at  opposite  sides  of  the  stack.  Under  a  simplified  assumption  for  the 
flow  of  oxygen  entering  the  gas  diffusion  layer  of  each  cell,  detailed  mass  flow  and  pressure  distributions 
are  shown,  including  the  possibility  of  a  turbulent  flow  inside  the  main  collectors. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

A  correct  distribution  of  the  reactant  flow  through  each  one 
of  the  cells  that  make  up  a  proton  exchange  membrane  fuel  cell 
(PEMFC)  stack  is  essential  to  a  correct  performance  of  the  fuel  cell. 
In  particular,  a  similar  distribution  of  mass  flow  circulating  through 
each  individual  cell  is  desirable,  which  also  helps  to  guarantee  a 
most  important  issue :  an  approximately  equal  pressure  distribution 
for  all  plates  [1  ].  It  should  be  reminded  that  the  velocity  of  the  gas 
reactants  entering  the  gas  diffusion  layer  (GDL)  is  essentially  depen¬ 
dent  on  the  pressure  gradient  along  it,  which  on  its  turn  depends  on 
the  geometry  and  flow  conditions.  As  all  the  plates  in  the  stack  have 
the  same  geometry,  then  it  is  the  flow  conditions  which,  under  an 
ideal  situation  (no  water  flooding,  etc.),  finally  determine  the  pres¬ 
sure  gradient,  and  eventually,  the  flow  of  reactants  reaching  the 
catalyst  layers. 

Some  previous  works  [2,3]  have  addressed  this  problem  by 
considering  a  standard  pipe  network  approach,  using  for  example 
friction  coefficients  to  calculate  both  the  pressure  losses  along  the 
main  collectors  and  local  losses  associated  to  bends.  The  main  crit¬ 
icism  for  this  otherwise  fast  approximation  is  that  the  “exits”  from 
the  collectors  to  the  individual  plates  are  multiple  and  very  close. 
That  disturbs  the  flow  in  the  main  collectors,  so  for  example,  in  a 
laminar  case  in  a  cylindrical  collector,  the  velocity  does  not  present 
the  typical  paraboloidal  Poiseuille  shape.  Something  alike  happens 
in  the  case  of  a  turbulent  flow  inside  the  collectors.  The  consequence 
is  that  the  real  friction  coefficient  does  not  match  the  theoretical 
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one.  In  fact,  the  real  coefficient  depends  strongly  on  the  geometry 
and  flow  conditions,  making  this  approach  dubious. 

Another  more  recent  approach  [4]  considers  a  2D  approximation 
where  the  channels  are  not  represented  geometrically,  but  taken  as 
tubes  filled  out  with  porous  media. 

In  this  paper,  the  Navier-Stokes  (NS)  equations  are  solved  in 
the  main  collectors  (with  a  turbulence  model)  and  channels  of  the 
plates,  with  a  detailed  3D  description  and  without  any  extra  model 
for  the  flow  or  geometry.  The  flow  through  the  porous  media  (GDL) 
is  also  solved  using  the  Brinkman-Darcy  approximation,  as  the  for¬ 
mulation  by  Ochoa-Tapia  and  Whitaker  [5]. 

2.  Mathematical  formulation 


As  explained  in  Section  1,  to  specify  flow  and  pressure  distribu¬ 
tions  along  the  stack,  including  the  plates,  NS  equations  are  solved, 
with  the  needed  modifications  for  the  flow  in  the  porous  media. 
In  the  case  of  a  laminar  flow,  NS  equations  are  fully  solved  and  no 
extra  model  is  used.  In  the  case  of  a  turbulent  flow,  the  Standard 
k-e psilon  model  is  used.  Steady-state  conditions  are  considered  in 
both  laminar  and  turbulent  configurations. 

In  the  channels  of  the  bipolar  plates  and  the  collectors,  the  3D 
steady  version  of  the  incompressible  Navier-Stokes  equations  is 
used,  v  being  the  kinematic  viscosity,  p  the  density  and  ut  the  i  = 
1,2,3  component  of  the  velocity  field: 

Continuity : 


Momentum : 

dui  1  dp  (  d1 2ut 
Uj  dxj  ~  p  dxt  +  V  dxj  dxj  ’ 


(1) 

(2) 
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For  the  GDL,  it  is  convenient  to  distinguish  between  two  types  of 
volume  averages  of  the  physical  magnitudes  of  interest  inside  the 
porous  media  [5]: 

Superficial  average : 

<3) 


Intrinsic  average: 


where  j3  indicates  the  zone  available  for  the  fluid  inside  the  porous 
media  (the  “void”  part  in  opposition  to  the  solid  part),  and  V  indi¬ 
cates  a  small  enough  volume  used  to  calculate  the  average.  In 
agreement  with  the  notation  used,  Vp  indicates  the  void  part  inside 
the  volume  V,  its  fraction  being  the  porosity  s  by  definition.  Hence 
superficial  and  intrinsic  averages  are  related  through  (•)  =  £{■)&. 

It  can  be  shown  [  5  ]  that  under  appropriate  conditions,  the  super¬ 
ficial  averaged  velocity  is  the  matching  quantity  to  the  flow  velocity 
and  that  the  intrinsic  averaged  pressure  is  the  matching  quantity  to 
the  flow  pressure  inside  the  porous  media.  However  to  simplify  the 
notation,  u  and  p  will  be  directly  written  for  equations  inside  the 
porous  media.  For  a  steady  state  approximation,  the  conservation 
equations  are 

Continuity  inside  GDL : 


Momentum  inside  GDL : 

1  3 ut  _  1  dp  v  d2Uj  v 

s2  Uj  dxj  ~  p  dxt  +  s  dxj  dxj  I<  U] 


(5) 

(6) 


Fig.  1.  Transversal  cut  through  the  mid-plane  of  a  bipolar  plate  of  the  fuel  cell  stack 
depicting  the  cascade  design  of  the  channels.  This  is  common  to  all  the  bipolar  plates 
that  were  used  for  simulations. 

2.1.  Computational  parameters 

2.1.1.  Numerics 

Eqs.  (5)  and  (6)  have  been  discretized  with  a  finite  volume 
method  on  a  4  M  cell  hexahedral  mesh.  The  numerical  calculations 
have  been  carried  out  using  OpenFOAM  software  (version  1.4.1), 
with  the  SIMPLE  scheme.  For  solving  the  linear  systems  of  equa¬ 
tions,  a  bi-conjugate  gradient  scheme  for  velocities  and  an  algebraic 
multi-grid  one  for  pressure  have  been  used.  The  code  is  ran  in  a 
parallel  environment  through  domain  decomposition.  The  conver¬ 
gence  time  on  a  Beowulf  with  8  Opteron  (64  bits  processor)  and 
96  GB  RAM  lasted  no  longer  than  a  few  hours. 


where  s  is  the  porosity  and  I<  is  the  permeability,  assuming  an 
isotropic  and  homogeneous  porous  medium.  The  second  term  on 
the  right  hand  side  of  Eq.  (6)  is  known  as  the  Brinkman  approxima¬ 
tion,  and  the  third  one  reflects  the  contribution  of  Darcy’s  law.  One 
can  observe  that  Eq.  (2)  are  recovered  in  the  case  of  porosity  in  the 
channels  equal  to  unity  and  infinite  permeability.  Using  simplifying 
assumptions  [5],  Eqs.  (5)  and  (6)  can  be  used  throughout  the  whole 
domain  without  the  need  of  any“internal”  boundary  conditions. 

In  this  paper  the  configuration  considered  consists  of  a  stack 
with  the  collectors  and  the  cathodic  side  of  24  bipolar  plates.  In 
order  to  solve  the  equations,  appropriate  boundary  conditions  are 
to  be  prescribed  at  the  inlet  and  outlet  of  the  collectors  as  well  as 
at  the  outlet  of  each  cathodic  GDL.  Actually,  the  flow  at  the  lat¬ 
ter  boundaries  is  affected  by  the  other  components  of  the  fuel  cell. 
However,  in  an  ideal  case,  the  stoichiometric  flow  can  be  considered 
as  a  good  approximation  and  this  can  be  imposed  by  specifying  a 
constant  (and  equal  in  magnitude  for  all  plates)  pressure  boundary 
condition  at  a  certain  distance  from  the  end  of  each  GDL.  Further¬ 
more,  in  the  case  studied,  only  oxygen  is  considered  as  a  reacting 
gas,  meaning  that  following  the  usual  practice  of  injecting  two  times 
the  stoichiometric  flow  and  the  fact  that  oxygen  is  only  a  fifth  of 
volume  air,  the  flow  at  the  outlet  collector  -  the  one  which  has 
not  reacted  -  would  be  one  order  of  magnitude  bigger  than  the 
sum  of  the  outlet  flows  at  the  GDL  exits— the  ones  consumed  by 
reaction.  Therefore  this  approximate  pressure  condition  at  the  GDL 
exits  should  not  greatly  affect  the  flow  through  the  collectors  and 
the  cathodic  channels.  Obviously,  the  simulation  assumes  a  stable 
operation  of  the  fuel  cell  stack. 

Notice  that  this  approximation  differs  from  the  ones  using  a 
pipe-network  approach,  as  the  flow  is  allowed  to  enter  the  GDL, 
which  is  much  more  realistic. 


2.1.2.  Design  and  physical  conditions 

As  mentioned  above,  the  parts  of  the  PEMFC  simulated  are  the 
collectors  and  the  cathodic  side  of  24  bipolar  plates.  Each  bipolar 
plate  is  designed  to  obtain  a  cascade-type  flow  pattern,  as  shown 
in  the  cut  throug  its  mid-plane  (see  Fig.  1),  with  an  active  area  of 
5  cmx  5  cm,  that  is  2500  mm2.  The  collectors  are  cylindrical  in 
shape  with  an  inner  diameter  of  6  mm  and  the  external  I/O  stack 
connectors,  also  included  in  the  domain,  have  an  inner  diameter  of 
4  mm.  The  connectors  in  betweeen  the  collectors  and  each  bipolar 
plate  (input  and  exit  ducts  for  each  plate)  have  a  cross  section  of 
1  mm2and  the  channel  height  in  the  bipolar  plates  is  of  1  mm.  The 
GDL  thickness  is  0.3  mm.  The  gas  is  air  at  80  C,  which  has  a  density 
p  =  0.996  kg  m-3  and  a  kinematic  viscosity  v  =  2.05  x  10-5  m2  s-1 . 
The  porosity  of  the  GDL  is  taken  to  be  e  =  0.75  with  a  permeability 
I<  =  8.69  x  10-12  m2.  The  flow  outlet  velocity  is  fixed  at  0.9 x  flow 
inlet  velocity,  following  the  reasonings  given  above. 

The  size  and  number  of  bipolar  plates  has  been  chosen  such  that 
the  cell  will  produce  12  V  considering  a  current  density  of  1  A  cm-2, 
each  plate  giving  0.5  V. 

The  cases  that  have  been  studied  are  synthesized  in  Table  1.  The 
total  power  of  the  stack  is  300  W  (Cases  II  and  III).  For  Case  I,  a  flow 


Table  1 

Numerically  simulated  cases.  Case  I  corresponds  to  a  laminar  flow  in  the  collectors 
while  Cases  II  and  III  to  turbulent  flows. 


Case 

Configuration 

Inlet  velocity  (m  s-1 ) 

Re 

I 

U 

3.4067 

664.7 

II 

U 

34.067 

6647 

III 

Z 

34.067 

6647 
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Fig.  2.  Comparisons  between  numerical  and  experimental  data  of  pressure  drop 
between  the  entry  and  exit  of  a  bipolar  plate  for  different  flows. 


10-fold  smaller  was  considered,  so  that  laminar  conditions  in  the 
collectors  are  attained  and  studied. 

3.  Results  and  discussion 

3.1  Pressure  drop  at  the  bipolar  plates 

First  of  all,  and  in  order  to  test  the  numerical  code,  a  set  of  com¬ 
parisons  with  simple  experimental  data  has  been  carried  out  for  the 
flow  circulating  through  the  channels  of  a  plate,  for  the  “cascade” 
configuration  studied  in  this  paper.  To  measure  the  pressure  drop  in 
the  plate,  two  branches  of  a  differential  inclined-tube  manometer 
are  connected  to  the  inlet  and  outlet  ducts.  For  both  experimen¬ 
tal  measurements  and  numerical  simulations  no  flow  through  the 
GDL  is  considered.  As  Fig.  2  depicts,  there  is  an  excellent  agreement 
between  the  pressure  drop  measured  and  simulated,  even  for  the 
low  values  of  drop  detected  in  the  cascade  configuration,  which 
proves  the  good  behaviour  of  this  geometrical  pattern  (and  of  the 
numerical  schemes). 

It  could  be  at  first  surprising  the  low  pressure  drop  through  the 
bipolar  plates  compared  with  the  one  taking  place  at  the  entry  (and 
exit)  ducts  to  the  collectors  (see  the  following  figures).  Naturally,  the 
reason  lies  on  the  verly  low  value  of  the  velocity  of  the  air  circulating 
through  the  channels. 

3.2.  U  and  Z  configuration  stack 

The  laminar  situation  would  be  the  ideal  one  due  to  the  lower 
pressure  drops.  Fig.  3  depicts  a  transversal  cut  through  the  stack  in  U 
configuration  in  the  laminar  case  (I).  One  can  observe  that  the  pres¬ 
sure  profile  does  not  remain  uniform  along  the  collector  transversal 
sections,  meaning  that  the  velocity  profile  does  not  preserve  its 
paraboloidal  shape.  As  a  consequence,  the  friction  coefficient  for  a 
laminar  flow  in  a  circular  section  duct  cannot  be  applied. 

As  previously  mentioned  it  is  desirable  that  the  flow  in  the 
collectors  be  laminar.  However,  our  condition  for  the  flow  to  be 
injected  in  the  stack  geometry  imposes  relatively  high  Reynolds 
(Re)  numbers.  Thus,  as  mentioned  above,  some  kind  of  turbulence 
modeling  is  needed  for  Cases  II  and  III  and  we  chose,  in  the  RANS 
framework,  the  Standard  k-epsilon  model,  as  it  has  been  mentioned 
above. 

Note  that  at  low  Re,  this  model  recovers  the  NS  equations,  as 
required,  because  as  it  can  be  seen,  the  flow  in  the  connecting 
ducts  between  the  bipolar  plates  and  the  collectors,  as  well  as  the 
flow  through  the  bipolar  plates  themselves  is  laminar.  This  can  be 


Fig.  3.  Transversal  cut  for  the  U  configuration  of  the  stack  in  the  laminar  case  (I).  It 
is  observed  that  the  pressure  contour  pattern  in  the  collectors  is  not  consistent  with 
a  laminar  flow  in  a  duct. 


observed  in  Fig.  4,  where  it  can  be  seen  that  the  flow,  although  tur¬ 
bulent  in  the  main  part  of  the  collectors,  transforms  rapidly  into 
laminar  in  the  connecting  ducts  and  the  bipolar  plates  where  it  has 
a  typical  paraboloidial  shape.  As  it  is  expected  in  a  porous  media, 
there  is  a  sudden  pressure  drop  across  the  GDL.  It  can  be  seen 
in  Fig.  5  the  rapid  decay  in  the  mass  flow  in  the  inflow  collector 
due  to  the  multiple  outlets  towards  the  bipolar  plates.  The  typical 
paraboloidal  shape  for  the  velocity  field  in  laminar  flows  is  however 
recovered  inside  the  small  entry  (and  exit)  ducts  to  each  plate. 

The  pressure  drop  takes  place  mostly  at  the  entry  (and  exit)  ducts 
to  the  bipolar  plates.  This  is  due  both  to  the  90°  bend  and  to  the 
relatively  high  velocity  and  very  small  cross  section  of  those  ducts 
(see  Figs.  4, 6  and  7).  The  flow  in  the  bipolar  plates  is  almost  uniform 
and  with  a  very  small  pressure  drop,  just  what  it  is  to  be  desired  for 
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Fig.  5.  Flow  at  inlet  collector-bipolar  plate  for  the  Z  configuration  in  the  turbulent 
case.  Velocity  details. 
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Fig.  6.  Transversal  cut  of  the  stack  for  turbulent  flow  both  for  the  U  configuration  Case  II  (a),  and  Z  configuration  Case  III  (b). 


Fig.  7.  3D  detail  of  collector-connector-bipolar  plate-GDL  for  the  Z  configuration  and 
turbulent  case  (III). 

a  correct  operation  of  the  stack.  This  is  typical  to  the  cascade  design 
[1  ]  and  this  can  be  also  seen  in  Fig.  6  where  both  configurations  are 
presented. 

In  Fig.  7  a  3D  detail  of  the  flow  for  the  stack  in  Case  III  is  presented. 
One  can  observe  the  typical  pressure  drop  in  the  GDL.  This  effect 


Fig.  8.  Mass  flow  distribution  along  the  bipolar  plates  inlets  for  both  U  and  Z  con¬ 
figurations. 


is  relatively  small  due  to  the  small  values  of  the  component  of  the 
velocity  orthogonal  to  the  interface  between  the  channels  and  the 
porous  media. 

In  Fig.  8  the  mass  flow  rate  through  each  one  of  the  bipolar  plates 
is  presented  for  both  U  and  Z  configuration  in  the  turbulent  case  (II 
and  III).  In  principle  the  Z  configurations  provides  a  more  uniform 
mass  flow  distribution.  However,  as  it  can  be  noticed  in  the  fig¬ 
ure,  the  effect  is  not  as  important  as  to  prevent  the  use  of  the  U 
configuration  when  required. 

4.  Conclusions 

In  the  present  paper  it  has  been  demonstrated  the  possibility  and 
interest  of  numerical  simulations  of  the  full  air  flow  inside  a  300  W 
PEM  fuel  cell  under  simplified  working  conditions.  The  flow  of  air 
has  been  studied  for  the  whole  stack,  including  the  GDL,  consider¬ 
ing  full  3D  geometry  and  full  NS  equations,  with  no  flow  modelling 
assumptions,  other  than  a  Standard  k-espsilon  closure  for  the  tur¬ 
bulence  which  exists  in  the  main  collectors.  The  flow  in  the  GDL 
has  been  also  included  in  the  simulations,  using  a  Brinkman-Darcy 
approximation  for  porous  media. 

It  has  been  shown  that  the  limitations  of  pipe-network  approxi¬ 
mations  make  them  inappropriate  for  the  use  to  calculate  gas  flows 
in  PEMFC  stacks.  The  main  reason  is  the  large  number  of  exits  to 
(and  from)  the  plates,  which  forces  the  flow  in  the  collectors  not  to 
follow  the  typical  pattern  inside  a  cylindrical  duct. 

On  the  contrary,  the  full  NS  approximation  proves  to  be  adequate 
while  not  expensive  computationally,  allowing  a  detailed  in-space 
description  of  the  flow  and  pressure  fields.  For  example  it  can  be 
seen  that  most  of  the  pressure  drops  take  place  at  the  small  ducts 
entering  and  exiting  each  plate,  not  in  the  GDL,  the  channels  of  the 
plates  or  the  main  collectors. 

Two  basic  configurations  with  turbulent  flow  in  the  collectors 
have  been  studied  (U  and  Z).  To  complete  the  study,  a  full  laminar 
configuration  in  U  has  been  also  studied.  In  all  cases  the  flow  pattern 
is  adequate,  although  the  Z  configuration  is  preferable. 
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